﻿using System;

namespace SoftwareConsulting.BI11.NumericalMethods.Optimizing
{
    /// <summary>
    /// Класс минимизации функции методом Брента
    /// </summary>
    public class BrentMethod
    {
        private const double _cGold = 0.3819660;
        private double _resultMinimumPoint; //точка найденного минимума
        private double _result; //минимум
        
        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="function">функция, требующая решения</param>
        /// <param name="leftBorder">левая граница отрезка, на котором ищется минимум</param>
        /// <param name="rightBorder">правая граница отрезка, на котором ищется минимум</param>
        /// <param name="epsilon">абсолютная точность, с которой находится расположение минимума</param>
        public BrentMethod(FunctionWithOneArgument function, double leftBorder, double rightBorder, double epsilon)
        {                                 
            double bx = 0.5 * (leftBorder + rightBorder);
            double ia = leftBorder < rightBorder ? leftBorder : rightBorder;
            double ib = leftBorder > rightBorder ? leftBorder : rightBorder;                        
            double v = bx;
            double w = v;
            double x = v;
            double e = 0.0;
            double fx = function(x);
            double fv = fx;
            double fw = fx;
            double d = 0;
            double u = 0;
            Func<double, double, double> signDefinder = (a, b) => b > 0 ? Math.Abs(a) : -Math.Abs(a);                
            for (int iter = 1; iter <= 100; iter++)
            {
                double xm = 0.5 * (ia + ib);
                if (Math.Abs(x - xm) <= epsilon * 2 - 0.5 * (ib - ia))
                    break;                
                if (Math.Abs(e) > epsilon)
                {
                    double r = (x - w) * (fx - fv);
                    double q = (x - v) * (fx - fw);
                    double p = (x - v) * q - (x - w) * r;
                    q = 2 * (q - r);
                    if (q > 0)
                        p = -p;                    
                    q = Math.Abs(q);
                    double etemp = e;
                    e = d;
                    if (!(Math.Abs(p) >= Math.Abs(0.5 * q * etemp) | p <= q * (ia - x) | p >= q * (ib - x)))
                    {
                        d = p / q;
                        u = x + d;
                        if (u - ia < epsilon * 2 | ib - u < epsilon * 2)
                            d = signDefinder(epsilon, xm - x);                        
                    }
                    else
                    {
                        e = x >= xm ? ia - x : ib - x;
                        d = _cGold * e;
                    }
                }
                else
                {
                    e = x >= xm ? ia - x : ib - x;                    
                    d = _cGold * e;
                }
                u = x + (Math.Abs(d) >= epsilon ? d : signDefinder(epsilon, d));                
                double fu = function(u);
                if (fu <= fx)
                {
                    if (u >= x)
                        ia = x;                    
                    else
                        ib = x;                    
                    v = w;
                    fv = fw;
                    w = x;
                    fw = fx;
                    x = u;
                    fx = fu;
                    continue;
                }                
                if (u < x)
                    ia = u;                
                else
                    ib = u;                
                if (fu <= fw | w == x)
                {
                    v = w;
                    fv = fw;
                    w = u;
                    fw = fu;
                    continue;
                }                
                if (fu <= fv | v == x | v == 2)
                {
                    v = u;
                    fv = fu;
                }                                
            }
            _resultMinimumPoint = x;
            _result = fx;
        }

        /// <summary>
        /// Возвращает решение - точку минимума
        /// </summary>        
        public double GetSolutionMinimumPoint()
        {
            return _resultMinimumPoint;
        }
        
        /// <summary>
        /// Возвращает решение - минимум функции на указанном отрезке
        /// </summary>
        /// <returns></returns>
        public double GetSolutionMinimum()
        {
            return _result;
        }
    }
}
